UNIVEESITY 


Center  for  Air  Sea  Technology 


Prepared  for:  Office  of  Naval  Research  . 

Under  Grants  N0001 4-97-1 -0525  ^  ^  UySPEGTOp  q 

and  N00014-97-1-0099 


Approved  for  public  release;  distribution  is  unlimited. 
Mississippi  State  University  Center  for  Air  Sea  Technology 
_ Stennis  Space  Center,  MS  39529-6000 _ 


TECHNICAL  REPORT  03-98 


A  SEMI-IMPLICIT  FREE  SURFACE  FORMULATION 
FOR  THE  SEMI-COLLOCATED  GRID 
DIECAST  OCEAN  MODEL 


by 


DAVID  E.  DEITRICHi  and  AVICHAL  MEHRA2 


1  Senior  Research  Scientist,  Mississippi  State  University  Center  for  Air  Sea  Technology,  Stennis 
Space  Center,  MS  39529-6000 

2  Postdoctoral  Assistant,  Mississippi  State  University  Center  for  Air  Sea  Technology,  Stennis 
Space  Center,  MS  39529-6000 


15  APRIL  1998 


This  research  was  supported  by  the  Department  of  the  Navy,  Office  of  Naval  Research  under 
Grants  NOOOl  4-97- 1-0525  and  NOOO 14-97- 1-0099  with  Mississippi  State  University.  The 
opinion,  findings,  conclusions,  and  recommendations  expressed  in  this  publication  are  those  of 
the  authors  and  do  not  necessarily  reflect  the  views  of  the  U.S.  Government.  No  official 
endorsement  should  be  inferred. 


Table  of  Contents 

Page  Number 

Cover  Sheet . i 

Table  of  Contents  . ii 

Abstract . iii 

1.1  INTRODUCTION .  1 

1.2  NEW  SEMI-IMPLICIT  FREE  SURFACE  FORMULATION 

FOR  THE  SEMI-COLLOCATED  DIECAST  MODEL . 2 

1.3  PROCEDURE  SUMMARY . 5 

1.4  COMMENTS  ON  PRESSURE  GRADIENT  EVALUATION  IN 

SEMI-IMPLICrr  MIXED  “a”  and  “c”  GRID  MODELING . 6 

1.5  FREE-SURFACE  RESULTS .  7 

ACKNOWLEDGEMENTS 

REFERENCES 
DISTRIBUTION  LIST 


REPORT  DOCUMENTATION 


A  SEMI  -  IMPLICIT  FREE  SURFACE  FORMULATION  FOR 
A  SEMI-COLLOCATED  GRID  DIECAST  OCEAN  MODEL 


ABSTRACT 

We  present  a  semi-implicit  free-surface  formulation  for  the  DieCAST 
Ocean  model  that  retains  the  accurate,  low  dissipation  numerics  of  its  lat¬ 
est  and  best  semi-collocated  rigid-lid  version.  The  approach  involves  inte¬ 
grating  the  time-implicit  (trapezoidal)  shallow  water  equations  on  a  stag¬ 
gered  Arakawa  “c”  grid,  with  vertically  averaged  baroclinic  forcing  terms 
determined  on  a  non-staggered  Arakawa  “a”  grid,  including  a  fourth-order- 
accurate  baroclinic  pressure  gradient.  The  shallow  water  equation  numer¬ 
ics  are  virtually  equivalent  to  standard  sigma  coordinate  approaches  for  the 
model  barotropic  mode.  In  application  to  transient  wind  forced  lake  Kelvin 
waves,  the  new  free-surface  version  gives  results  very  close  to  the  correspond¬ 
ing  strongly  validated  rigid-lid  DieCAST  version  (reflecting  the  fact  that  the 
rigid-lid  barotropic  mode  numerics  are  also  a  sigma-like  approach),  and  re¬ 
quires  less  than  50  percent  more  computing.  Thus,  the  numerics  used  by  both 
rigid-lid  and  free-surface  DieCAST  versions  combines  the  best  of  z-level  and 
sigma  coordinate  numerics  as  well  as  the  best  of  “a”  and  “c”  grid  numerics. 


1.1  Introduction 

In  many  computational  fluid  dynamics  applications,  the  flows  are 
steady  or  slowly  changing  compared  to  the  shortest  time  scale  resolved  modes. 
In  such  case,  it  may  be  desirable  to  take  longer  time  steps  than  explicit  meth¬ 
ods  (e.g.,  leap  frog)  allow,  especially  when  nonlinear  low  frequency  effects  (if 
any)  of  the  high  frequency  modes  are  parameterized  (i.e.,  given  by  ad  hoc 
relations  to  the  lower  frequency  modes).  Such  parameterization  may  also  in¬ 
clude  the  effects  of  subgrid-scale  modes  that  are  not  resolved  either  spatially 
or  temporaly. 

Semi-implicit  and  fully-implicit  time  marching  methods  allow  these 
longer  time  steps  (e.g.,  Kwizak  and  Robert,  1971;  O’Brien  and  Hurlburt, 
1972;  Dietrich,  et.  al,  1975;  Dietrich,  1975;  Dietrich  and  Wormeck,  1985). 
These  methods  generally  “slow  down”  the  fastest  (highest  frequency)  .modes, 
but,  as  noted  above,  this  may  be  acceptable.  Indeed,  the  fastest  free  surface 
modes  have  only  secondary  effect  on  the  ocean  general  circulation.  This  was 
necessary  for  the  success  of  popular  ocean  models  that  use  a  rigid-lid  approx¬ 
imation  or  semi-implicit  calculation  of  the  high  frequency  surface  modes.  In 
converting  the  popular  rigid-lid  Bryan-Cox  model  to  a  corresponding  model 
having  semi-  implicit  free  surface  waves,  Dukowicz,  et.  al  (1992)  showed  that, 
after  initial  transients,  semi-implicit  free  surface  models  ocean  general  circu¬ 
lation  results  are  virtually  identical  to  rigid-lid  models,  rather  independent 
of  the  time  step  size  up  to  the  rigid- lid  stability  limit. 

Semi-implicit  and  rigid-lid  models  generally  require  solution  of  an 
elliptic  equation  each  time  step.  Thus,  a  good  elliptic  solver  is  needed  in 
order  for  them  to  be  useful.  A  major  advantage  of  hydrostatic  ocean  models 
is  that  the  hydrostatic  approximation  reduces  the  elliptic  equation  to  two 
dimensions,  thereby  greatly  reducing  elliptic  solver  requirements  (compared 
to  non-hydrostatic  three-dimensional  models)  and  allowing  use  of  marching 
“error  vector  propagation”  (EVP)  methods  (Roache,  1995)  that  are  highly 
efficient  for  low-  and  moderate-  resolution  grids. 

However,  for  very  large  grids  (greater  than  0(500)  points  in  the  “cross¬ 
march”  direction)  EVP  storage  requirements  can  be  large  and  sometimes 
EVP  cannot  take  full  advantage  of  vector  and  massively  parallel  supercom¬ 
puters  (vectorization  and  parallelization  is  usually  possible  only  in  the  cross¬ 
march  direction).  In  such  case,  it  may  be  useful  to  combine  EVP  with  do¬ 
main  decomposition  (Dietrich,  1974;  Dietrich,  et.  al,  1975;  Roache,  1995), 
thus  leading  to  an  iterative  solution  procedure.  Iterative  procedures  work 
best  for  diagonally  dominant  elliptic  problems.  The  elliptic  problem  asso¬ 
ciated  with  rigid-lid  ocean  models  is  a  Poisson  equation  (Dietrich,  et.  al, 
1987;  Dietrich,  1992)  which  is  less  diagonally  dominant  than  the  Helmholtz 
equation  associated  with  semi-implicit  free  surface  models. 
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It  follows  that,  besides  being  able  to  address  free-surface  modes,  semi- 
implicit  free  surface  ocean  model  formulations  may  be  useful  numerically  for 
very  large  grids  even  when  the  free  surface  modes  are  of  secondary  interest. 
Thus,  although  low-to- moderate  resolution  rigid-lid  models  may  remain  at¬ 
tractive  for  general  circulation  simulations,  high  resolution  grids  may  favor 
semi-implicit  free  surface  formulations. 

Previously,  the  DieCAST  ocean  model  has  used  a  rigid-lid  approxima¬ 
tion.  Now  that  the  new  rigid-lid  semi-collocated  versions  have  been  validated 
by  extensive  Gulf  of  Mexico  observations  (Dietrich  and  Lin,  1997;  Dietrich, 
1997)  we  are  herein  led  to  formulate  a  semi-implicit  free  surface  version  that 
retains  the  demonstrated  accurate,  robust,  low  dissipation,  low  dispersion 
numerics  of  the  latest  and  best  rigid-lid  versions. 

1.2  New  Semi-Implicit  Free  Surface  Formulation  for  the  Semi- 
Collocated  DieCAST  Ocean  Model 

A  major  reason  for  the  success  of  the  DieCAST  ocean  model  has  been 
its  robust,  accurate,  low  dispersion  and  low  dissipation  numerics.  We  now 
formulate  a  semi-implicit  free  surface  version  that  includes  its  latest  and  best 
numerics  (Dietrich,  1997). 

First,  it  is  convenient  to  write  the  governing  time-implicit  (Crank 
Nicolson  or  trapezoidal  approximation)  equations  in  terms  of  the  two-  level 
time  averaged  fields: 

=  jC?'  +  9’-"^')  (1) 

In  the  following  equations,  all  time-explicit  terms  are  given  in  their 
analytic  form.  They  are  calculated  using  the  standard  DieCAST  numerics, 
including  filtered  leapfrog  time  integration  and  fourth  order  advection  as 
before.  The  leapfrog  time  level,  at  which  all  explicit  terms  are  evaluated,  is 
at  time  ^  -I- 

The  incompressibility  (mass  conservation)  equation  is  solved  only  on 
an  Arakawa  “c”  grid  (see  below).  The  horizontal  momentum  conservation 
equations  to  be  solved  on  an  Arakawa  “a”  grid  are: 

~  “  9^  ”  ^ 

V^  =  V*+  ■  (vV  -  uVv)  -  -^-fu  +  Ty]  (3) 

where  the  wind  stress  Tj,,  Ty  applies  to  top  layer  only.  There  are  also  transport 
equations  for  temperature  and  salinity,  but  these  are  not  involved  in  the  semi- 
implicit  algorithm  described  here. 
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The  solution  procedure  is  as  follows.  First,  we  partially  update  (with¬ 
out  surface  height  or  baroclinic  pressure  gradient  terms)  the  semi-collocated 
(Arakawa  “a”  grid)  fields  u*,  u*  as  follows: 

u  =  -f  ^[- V  •  {uV  -  vVu)  +  fv  +  T^]  (4) 

V  =  V*  +  ^[—V  •  {vV  -  uVv)  —  fu  +  Ty]  (5) 

where  ii  =  +  f  (jf  +  g)  and  6  =  S*  +  f  (jf  +  g). 

Next,  we  interpolate  (to  4th  order  accuracy)  u  and  v  to  the  “c”  grid 
locations: 

—  interpolated  (u)  and  v‘  =  interpolated  (u)  (6) 

Next,  we  add  the  baroclinic  pressure  gradient  (standard  second  order 
accurate  approximation  is  adequate  on  the  “c”  grid): 

^  At  dp  At  dp  .  . 


where  p  =  Jq  pgdz  is  the  hydrostatic  baroclinic  pressure  head  based  on  the 
known  density  anomaly  (deviation  from  horizontal  average),  p  . 

Equations  (7)  include  all  time-explicit  terms  in  the  implicit  equations 
for  “c”  grid  velocity.  Adding  the  implicit  surface  height  gradient  terms  to 
the  right  hand  sides  of  equations  (7)  will  provide  the  fully  updated  “c”  grid 
two-level-time-averaged  velocity.  We  determine  the  surface  height  implicitly 
as  follows. 


Defining, 
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dz  and  V 


-r-‘ 

H  Jo 


where  H  is  the  local  depth  ignoring  small  surface  wave  fluctuations,  the 
implicit  “c”  grid  barotropic  mode  spatially  differenced  equations  are: 

+  Ui,  (9) 


3 


I  _  tt  r  ^i+l/2,j  ^i,j  ~  12, j  ,  Hi,j^il2  Vi,j  Hij-i/2  fij-i  ^ 

2  ^  Aa:  Ay  ^ 

(11) 

where  u  and  v  are  vertically  and  two- level  time  averaged  “c”  grid  quantities, 
not  to  be  confused  with  u*  and  in  the  three-dimensional  equations  (2)  and 
(3).  Equations  (9)-(ll)  are  an  implicit  time  differenced  formulation  of  the 
shallow  water  equations.  These  are  fundamentally  an  application  of  sigma 
coordinates  for  the  barotropic  mode. 

Substituting  equations  (9)  and  (10)  into  the  right  hand  side  of  equa¬ 
tion  (11)  gives  us  a  Helmholtz  equation  for  h  with  sources  from  U  and  V: 

^t-H/2,j(fet+l,j  ~  ^»',j)  ~  ^;'-l/2,j(^tJ  ~  ^1-1, i) 

(Aa;)2 

(Ay)2 


y(At) 


2  ~ 


^  Lt  ,  ^  /"^»+l/2,i  ^i,j  Hi-if2,j  Ui-i,j 

+  TTli- - 


y(Ai)2  gAV  Ax 

Hi,j+i/2  Vi,j  —  Hi,j_i/2  Vjj-i 


Ay 


-) 


(12) 


Equation  (12)  is  solved  for  h.  This  result  may  be  substituted  into 
equations  (9)  and  (10)  to  get  barotropic  “c”  grid  u  and  v,  but  these  are  not 
really  needed  in  this  formulation. 

Next,  the  full  3-d  “c”  grid  advanced  time  level  fields  are  updated  from: 


'>3“  =  2  Sy  -  Ki 


hij) 


=  2<,,,  -  sAf 


—  V: 


w 


Ay 

-  ■  Jh^ — a; — + — A^ — ^ 


(13) 

(14) 

(15) 

(16) 


( 
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Equations  (14-16)  give  a  non-divergent  3-d  advection  velocity  for  the 
next  time  step.  The  vertical  velocity  does  not  vanish  at  the  top  as  it  does 
in  rigid-lid  formulations.  Neuman  conditions  may  be  used  to  get  reasonable 
advective  fluxes  through  the  open  top,  because  the  field  variables  generally 
are  nearly  constant  over  depths  much  larger  than  h  variations.  A  more  accu¬ 
rate  approach  would  include  time  dependent  control  volumes  in  the  top  layer, 
but  the  present  simpler  approach  should  be  adequate  except  in  shallow  water 
and  large  amplitude  wave  conditions  (i.e.,  h  not  small  compared  to  top  layer 
thickness).  Large  dissipation  occurs  under  such  shallow  water  conditions,  al¬ 
lowing  appropriate  parameterization  through  model  dissipation  parameters; 
alternatively,  a  more  specialized  surf  zone  model  could  be  coupled  to  the 
present  model. 

Finally,  the  4th  order  “a”  grid  pressure  and  surface  height  gradients 
are  calculated  and  substituted  into  equations  (2)  and  (3)  to  finalize  the  3-d 
“a”  grid  solution. 

The  long  term  averaged  flow  resulting  from  the  above  procedure  will 
have  extremely  small  barotropic  mode  divergence,  because  there  is  no  long 
term  trend  of  the  ocean  surface  height.  As  noted  above,  the  short  term 
fluctuations  (the  so-called  “fast  modes” )  have  little  effect  on  the  ocean  general 
circulation.  It  follows  that  the  above  solution  procedure  will  yield  virtually 
the  same  ocean  general  circulation  as  the  earlier  rigid-lid  versions  of  the 
DieCAST  ocean  model  and  thus  retain  its  robust,  accurate,  low  dissipation 
and  low  dispersion  numerics. 

1.3  Procedure  Summary 

A  semi-implicit  free-surface  procedure  that  includes  the  latest  and 
best  4th  order  numerics  in  the  semi-collocated  DieCAST  version  has  been 
given.  The  procedure  is  as  follows: 

1.  determine  advection,  diffusion  and  Coriolis  terms,  and  wind  stress 
on  the  “a”  grid 

2.  interpolate  the  results  from  step  1  to  the  “c”  grid 

3.  evaluate  the  baroclinic  pressure  gradient  on  the  “c”  grid  using  the 
standard  second  order  accurate  “c”  grid  scheme  and  add  it  to  the  interpolated 
terms  from  the  “a”  grid 

4.  vertically  average  the  step  3  result 

5.  solve  the  implicit  trapezoidal  (Crank  Nicolson)  shallow  water  equa¬ 
tions  forced  by  the  step  4  result  (combined  to  form  a  Helmholtz  equation  for 
the  free  surface  height  that  is  identical  to  the  rigid  lid  pressure  equation  in 
the  limit  of  large  time  step) 

6.  update  the  three-dimensional  “c”  grid  velocity  by  adding  the  free- 
surface  height  gradient  calculated  in  step  5  to  the  terms  calculated  in  step  3. 
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7.  evaluate  the  full  pressure  gradient  on  the  “a”  grid  using  the  updated 
free  surface  height  calculated  in  step  5  and  adding  the  baroclinic  pressure 
gradient  (evaluated  to  4th  order  accuracy  on  the  “a”  grid)  (this  is  probably 
better  than  interpolating  “c”  grid  results  back  to  the  “a”  grid) 

8.  finish  the  time  step  using  results  of  steps  1  and  7. 

1.4  Comments  on  Pressure  Gradient  Evaluation  in  Semi-implicit 
Mixed  “a”  and  “c”  Grid  Modeling 

Two  possible  procedures  regarding  the  pressure  gradient  evaluation 
in  semi-implicit  free-surface  mixed  “a”  and  “c”  grid  models  are: 

1.  the  above  procedure,  which  evaluates  both  the  baroclinic  and  the 
barotropic  pressure  gradient  to  second-order  accuracy  on  the  natural  stag¬ 
gered  “c”  grid  of  the  barotropic  mode  equations,  and  also  evaluates  both  to 
fourth-order  accuracy  in  the  collocated  “a”  grid  equations; 

2.  a  modified  procedure,  which  evaluates  only  the  barotropic  (sea 
surface  height)  pressure  gradient  terms  on  the  “c”  grid  (to  second-order  ac¬ 
curacy)  in  the  implicit  barotropic  mode  equations.  Instead  evaluating  the 
baroclinic  pressure  gradient  on  the  “c”  grid,  the  fourth-order  baroclinic  “a” 
grid  pressure  gradient  is  interpolated  to  the  “c”  grid  (again  using  fourth-order 
accurate  interpolations). 

While  the  first  procedure  is  aesthetically  pleasing,  we  have  found  that 
it  is,  disappointingly,  unstable  when  using  long  implicit  time  steps  (more  than 
a  few  times  the  explicit  limit).  On  the  other  hand,  the  second  procedure  when 
carefully  done  to  full  fourth-order  accuracy,  is  more  accurate  than  the  first 
procedure  and  is  found  to  be  stable  with  very  long  time  steps,  not  limited 
by  explicit  limit.  For  example,  we  have  run  40  model  days  (1440  time  steps) 
with  time  step  size  ~50  times  the  explicit  limit. 

A  third  possible  procedure  is  to  evaluate  the  baroclinic  pressure  gra¬ 
dient  to  fourth-order  accuracy  on  the  “c”  grid  (together  with  a  second-order 
accurate  barotropic  pressure  gradient),  and  to  interpolate  the  mixed-order 
full  pressure  gradient  to  the  “a”  grid.  These  choices  apply  to  rigid-lid  as  well 
as  to  free-surface  approaches. 

The  barotropic  mode  pressure  (surface  height)  gradient  has  only  second- 
order  accuracy  in  all  three  procedures.  We  have  not  tested  the  third  proce¬ 
dure,  because  we  have  found  the  second  procedure  to  be  robust  as  well  as 
reasonably  accurate,  as  demonstrated  below  in  Section  1.5.  It  gives  accurate 
fast  mode  (barotropic  mode  gravity  wave)  phase  speed  when  the  fast  mode 
time  scale  is  resolved,  and  has  fourth-order  accuracy  on  the  baroclinic  modes 
that  dominate  ocean  general  circulation.  Although  the  second  procedure 
appears  best  at  this  time,  the  first  and  third  procedure  may  be  useful  when 
very  accurate  fast  mode  results  are  desired,  especially  in  ultra-high  resolution 
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coastal  simulations. 

Having  tested  the  second  procedure  in  idealized  applications,  we  are 
now  proceeding  toward  real  applications  and  iterative  elliptic  solver  imple¬ 
mentation  as  noted  above. 

1.5  Free-Surface  Results 

We  now  discuss  application  of  the  new  free-surface  DieCAST  version, 
and  the  reduced  dispersion  rigid-lid  DieCAST  “a”  grid  version  (Dietrich, 
1997)  to  a  transient- wind-forced  lake  Kelvin  wave  problem.  The  rigid- lid 
version  is  further  enhanced  by  suggestions  by  Dr.  Brian  Sanderson  (Univer¬ 
sity  of  New  South  Wales,  Australia)  and  Dr.  Dan  Wright  (Bedford  Institute 
of  Oceanography,  Canada). 

Specifically,  the  free-surface  and  rigid-lid  versions  are  both  applied 
with  625  m  resolution  to  an  idealized  circular  flat-bottom  lake  patterned 
after  summertime  Great  Lakes  conditions,  as  described  by  Beletsky,  et.  al 
(1997).  Both  versions  use  the  same  physical  parameters  and  wind  forcing, 
including  1  m^/sec  lateral  viscosity  and  diffusivity. 

Figure  1  shows  the  temperature  at  level  9.7  m  (near  the  thermocline) 
from  both  free-surface  and  rigid-lid  versions  at  days  6,  10,  14,  22  and  30. 
Although  the  results  differ  by  as  much  as  five  percent  during  the  first  few 
days,  they  agree  within  a  fraction  of  one  percent  after  day  ten.  This  is 
because  fast  free-surface  modes  disperse,  reflect  and  dissipate  during  the 
first  ten  days,  during  which  they  affect  thermocline  level  temperatures  as 
they  slosh  around  the  basin. 

Time  averaging  free-surface  results  over  1/2  day  to  reduce  fast  mode 
amplitudes,  and,  for  consistency,  applying  the  same  time  average  to  the  rigid- 
lid  results,  shows  in  Figure  2  that  the  early  results  are  closer  during  the  first 
few  days. 

Although  highly  inertial  Poincare  waves  are  strong  through  day  30, 
the  solution  is  dominated  by  the  Kelvin  edge  wave,  which  propagates  around 
the  100-km  diameter  lake  in  ~8  days. 

These  new  free-surface  and  rigid-lid  results  are  similar  to  results  from 
the  more  dispersive  original  “a”  grid  DieCAST  model  used  by  Beletsky  et. 
al  (1997).  However,  the  Kelvin  wave  speed  is  closer  to  the  theoretical  value 
and  the  results  are  more  strongly  inertial,  including  intense  internal  wave 
bore  development  after  day  20.  The  improvement  in  the  results  is  due  to 
accurate  reduced-dispersion  numerics  (Dietrich,  1997).  The  new  free-surface 
and  rigid-lid  DieCAST  results  are  much  closer  to  each  other  than  to  the 
original  “a”  grid  DieCAST  rigid-lid  results,  although  they  are  quite  similar 
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to  the  original. 

In  implicit  models,  increasing  time  step  size  reduces  phase  speeds  of 
modes  having  unresolved  frequencies.  However,  when  the  time  step  resolves 
the  wave  propagation,  the  new  semi-implicit  free-surface  DieCAST  version 
gives  accurate  free-surface  wave  phase  speeds  of  ~31  m/sec,  which  is  the 
correct  linear  phase  speed  for  the  100  m  deep  model  lake.  Its  fourth-order 
accuracy  on  the  baroclinic  terms  and  “c”  grid  accuracy  on  the  barotropic 
terms  guarantees  good  accuracy  when  resolution  is  adequate.  However,  a 
notable  advantage  is  that  the  time  step  can  be  many  times  longer  than  re¬ 
quired  to  resolve  the  highest  frequency  resolved  free  surface  “fast”  modes. 
This  does  not  destroy  the  accuracy  of  longer  wavelength  (and  thus  lower  fre¬ 
quency)  resolved  free  surface  modes.  Such  long  time  step  has  little  effect  on 
the  accuracy  of  slow  internal  waves  and  even  slower  baroclinic  modes,  or  on 
how  they  affect  free  surface  height. 

With  the  162  X  162  X  12  layers  resolution  used,  the  free-surface  ver¬ 
sion  requires  ~4.5  sec  cpu  per  time  step  on  an  SGI  Indigo  2  computer,  which 
runs  at  about  12,5  megaflops.  This  is  ~4.7  percent  more  computation  per 
time  step  than  the  rigid-lid  version.  The  rigid-lid  version  is  stable  with  time 
step  up  to  ~600  sec.  The  free-surface  version  is  stable  with  a  400  sec  time 
step,  and  unstable  with  a  600  sec  time  step.  With  a  400  sec  time  step,  the 
fast  surface  wave  Courant  number  is  ~50,  indicating  up  to  50  times  less  com¬ 
puting  compared  to  explicit  methods.  In  the  comparison  runs,  both  versions 
were  run  with  a  400  sec  time  step. 

It  is  noteworthy  that,  as  resolution  increases,  the  elliptic  solver  com¬ 
putation  and  storage  increase  relative  to  the  rest  of  the  model.  Thus,  for 
large  grids,  it  may  be  useful  to  explore  ways  to  reduce  the  elliptic  solver 
overhead.  As  noted  above,  a  significant  motivation  of  introducing  the  free 
surface  is  to  reduce  elliptic  solver  overhead  on  extra  large  grids,  because 
increased  diagonal  dominance  associated  with  the  free-surface  formulation 
benefits  potentially  lower  overhead  iterative  solvers  (Roache,  1995;  Dietrich, 
et.  al,  1975).  For  sufficiently  large  grids,  the  elliptic  solver  storage  over¬ 
head  is  substantially  reduced  and  its  computation  reduction  should  more 
than  compensate  for  the  slightly  smaller  time  step  required  by  the  DieCAST 
free-surface  version.  We  expect  the  free  surface  version,  combined  with  an 
iterative  elliptic  solver,  will  prove  especially  attractive  for  grids  larger  than 
0(500  X  500)  horizontal  resolution. 
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Figure  1:  Comparison  of  temperature  at  level  9.7  m  (near  the  thermocline) 
from  rigid-lid  (left)  and  free-surface  (right)  versions  at  days  6,  10,  14,  22  and 
30. 
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RisnlB(ire625m:  lemperature  (deg  c)  at  day  I0.dep0i-9.7m.  min- 1 1.83.  max-  13.69 


Figure  1:  (Coni.) 


HijulBcrcft’Sm:  Mnipmnir((de;e)aiday22.dep(M.7in.  min- 1 1.93.  max-  I3.0X 


FrireBorefi25m;  lempenture (degc) alday  22. depth- 9.7 m.  min-  1 1.94.  max= 


RigidBore625m:  l/2-<by-averaged  layer  4  lempentuer  at  day  6.  depth-  9.7  m.  min=  11.65.  max- 13.68 


Figure  2:  Comparison  of  time  averaged  results  (over  1/2  day)  from  rigid-lid 
(top)  and  free-surface  (bottom)  models  at  Day  6. 
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